
#include <fstream>
#include <string.h>
#include <cstdlib>
#include <iostream>
#include <vector>
#include <cmath>
#include "NewtonInterp.H"

const Real funB(const Real &_x)
{
    return 1.0 / (1 + _x * _x);
}

const Real funC(const Real &_x)
{
    return 1.0 / (1 + 25.0 * _x * _x);
}

const RealVect ChebyshevZeros(const int &N)
{
    RealVect temp;
    for (int i = 0; i < N; i++)
    {
        temp.push_back(cos((2 * i + 1) * M_PI / (2 * N)));
    }
    return temp;
}


using namespace std;
int main(int argc, const char **argv)
{
    cout << "-----------------------Assignment B-----------------------" << endl;
    NewtonInterp Solver;
    vector<Polynomial<Real>> polyBox;
    for (int n = 0; n < 4; n++)
    {
        int N = 2 * (n + 1);
        RealBox tempBox;
        InterpConditions IntpCondsB;
        for (int i = 0; i <= N; i++)
        {
            tempBox.push_back({-5 + 10.0 * i / N, funB(-5 + 10.0 * i / N)});
        }
        IntpCondsB.ResetIntpInfo(tempBox);
        Solver.ResetIntpConds(IntpCondsB);
        Solver.QuickInterp();
        polyBox.push_back(Solver.getInterPoly());
    }
    for (int n = 0; n < 4; n++)
        cout << polyBox[n] << endl;
    string bufB = "AssignmentB.m";
    char file_nameB[20];
    strcpy(file_nameB, bufB.c_str());
    ofstream fout(file_nameB);
    fout << "%This is a file about plotting AssignmentB." << endl;
    fout << "figure(1);" << endl;
    fout << "clear;x=-5:0.01:5;f=1./(1+x.*x);plot(x,f,'LineWidth',3);hold on" << endl;
    for (int n = 0; n < 4; n++)
    {
        fout << polyBox[n] << ";" << endl;
        fout << "plot(x,px,':','LineWidth',2);" << endl;
    }
    fout << "xlabel('x'),ylabel('y'),title('Runge Phenomenon');saveas(1,'AssignmentB.png');";
    fout.close();
    cout << "Generate AssignmentB.m!" << endl;
    polyBox.clear();

    cout << "-----------------------Assignment C-----------------------" << endl;
    for (int n = 0; n < 4; n++)
    {
        int N = 5 * (n + 1);
        RealBox tempBox;
        InterpConditions IntpCondsC;
        for (int i = 0; i < N; i++)
        {
            tempBox.push_back({ChebyshevZeros(N)[i], funC(ChebyshevZeros(N)[i])});
        }
        IntpCondsC.ResetIntpInfo(tempBox);
        Solver.ResetIntpConds(IntpCondsC);
        Solver.QuickInterp();
        polyBox.push_back(Solver.getInterPoly());
    }
    for (int n = 0; n < 4; n++)
        cout << polyBox[n] << endl;
    string bufC = "AssignmentC.m";
    char file_nameC[20];
    strcpy(file_nameC, bufC.c_str());
    ofstream fout1(file_nameC);
    fout1 << "%This is a file about plotting AssignmentC." << endl;
    fout1 << "figure(2);" << endl;
    fout1 << "clear;x=-1:0.01:1;f=1./(1+25*x.*x);plot(x,f,'LineWidth',3);hold on" << endl;
    for (int n = 0; n < 4; n++)
    {
        fout1 << polyBox[n] << ";" << endl;
        fout1 << "plot(x,px,':','LineWidth',2);" << endl;
    }
    fout1 << "xlabel('x'),ylabel('y'),title('Chebyshev Interpolation');saveas(2,'AssignmentC.png');";
    fout1.close();
    cout << "Generate AssignmentC.m!" << endl;
    polyBox.clear();

    cout << "-----------------------Assignment D-----------------------" << endl;
    InterpConditions IntpCondsD({{0,0,75},{3,225,77},{5,383,80},{8,623,74},{13,993,72}});
    Solver.ResetIntpConds(IntpCondsD);
    Solver.QuickInterp();
    polyBox.push_back(Solver.getInterPoly());
    cout<<polyBox[0]<<endl;
    polyBox.push_back(polyBox[0].derivate());
    cout<<"The derivation of "<<polyBox[1]<<endl;
    cout<<"(a) When t = 10, "<<endl;
    cout<<"Position: "<<polyBox[0].eval(10.0)<< " feet from Time 0."<<endl;
    cout<<"Speed: "<<polyBox[1].eval(10.0)<< " feet per second."<<endl;
    cout<<"(b) For t from 0 to 13, "<<endl;
    bool is_overspeed = false;
    for(int i=0; i<1300; i++)
    {
        if(polyBox[1].eval(0.01*i)>81.01)
        {
            is_overspeed = true;
            cout << "the car exceeds the speed limit when t = "<<0.01*i<<" ."<<endl;
            break;
        };
    }
    if(is_overspeed == false)
        cout << "the car never exceeds the speed limit. Innocent! "<<endl;
    polyBox.clear();

    cout << "-----------------------Assignment E-----------------------" << endl;
    InterpConditions IntpCondsE1({{0,6.67},{6,17.3},{10,42.7},{13,37.3},{17,30.1},{20,29.3},{28,28.7}});
    InterpConditions IntpCondsE2({{0,6.67},{6,16.1},{10,18.9},{13,15.0},{17,10.6},{20,9.44},{28,8.89}});
    Solver.ResetIntpConds(IntpCondsE1);
    Solver.QuickInterp();
    polyBox.push_back(Solver.getInterPoly());
    Solver.ResetIntpConds(IntpCondsE2);
    Solver.QuickInterp();
    polyBox.push_back(Solver.getInterPoly());
    for (int n = 0; n < 2; n++)
        cout << polyBox[n] << endl;
    string bufE = "AssignmentE.m";
    char file_nameE[20];
    strcpy(file_nameE, bufE.c_str());
    ofstream fout2(file_nameE);
    fout2 << "%This is a file about plotting AssignmentE." << endl;
    fout2 << "figure(3);" << endl;
    fout2 << "clear;x=0:0.01:28;hold on;" << endl;
    fout2 << "A"<< polyBox[0] << ";" << endl;
    fout2 << "plot(x,Apx,'LineWidth',2);" << endl;
    fout2 << "B"<< polyBox[1] << ";" << endl;
    fout2 << "plot(x,Bpx,'LineWidth',2);" << endl;    
    fout2 << "legend('Sp1','Sp2'),xlabel('x'),ylabel('y'),title('Average Weight Curves (0,28)');saveas(3,'AssignmentE1.png');"<<endl;
    fout2 << "figure(4);" << endl;    
    fout2 << "clear;x=28:0.01:43;hold on;" << endl;
    fout2 << "A"<< polyBox[0] << ";" << endl;
    fout2 << "plot(x,Apx,'LineWidth',2);" << endl;
    fout2 << "B"<< polyBox[1] << ";" << endl;
    fout2 << "plot(x,Bpx,'LineWidth',2);" << endl;    
    fout2 << "legend('Sp1','Sp2'),xlabel('x'),ylabel('y'),title('Average Weight Curves (28,43)');saveas(4,'AssignmentE2.png');";
    fout2.close();
    cout << "Generate AssignmentE.m !" << endl;
    polyBox.clear();
    return 0;
}


